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1 Introduction 

In a recent paper an extended Liouville- Green formula 

y{x)=jM{x){l + eM{x)}exp{± I QM[t)dt) (1.1) 



IX 

was developed for solutions of the second-order differential equation 

y"{x) - Q{x)y{x) = (0 < a; < oo). (1. 2) 

Here ^m{x) ~ Q~^{x), Qm(x) ^ Q^{x) and eM{x) — > as x ^ oo, while M(> 2) is an integer an d 7m 
and Qm can be defined in terms of Q and its derivatives up to order Af — 1. The general form of ( |l . l[ ) 
had been obta ined previously by Cassell (|] |^ and Eastha m pX | |pl] , section 2.4]. In particular, the 
proof of ( |l. l|) in ||l^ and |l^ depended on the formulation of ( 1. 2 ) as a first-order system and then on 



a process of M repeated diagonalizations of the coefficient matrices in a sequence of related differential 

systems. 

The main contribution of |^ to ( |1 . 1[ ) was to show that, for a general class of coefficients Q, the 
magnitude of the error term eM{x) for large x decreases as M increases. This feature of €m{x) was then 



exploited in ||3|] for the case where Q{x) = q{x) — A and (1. 2) becomes the usual Sturm-Liouville equation 
with spectral parameter A . The smallness of enix) for a suitable choice of M (such as M = 6) leads to 
an efficient numerical algorithm for estimating the Titchmarsh-Weyl function m(A) with precise global 
error bounds. In Q, the case where q{x) — — (0 < a < 2) and ReA > —1 was considered in detail. 

In this paper, we develop these ideas for higher-order differential equations. Thus we consider the 
generalisation of ( [L. 1[ ) for the n-th order equation 

y'-^Hx)-Q{x)y{x)^() (1.3) 



again with improving estimates on cm as M increases. In the case where n is even, n — 2v, and 
Q{x) = (— 1)'^{A — q{x)}, we have 'au xu spectral matrix (my (A)) which corresponds to the Titchmarsh- 
Weyl function m(A) p^ , and we can again apply our estimates for cm to obtain numerical estimates for 
the spectral matrix. 

There are however certa in difficulties which arise when n > 2 and which were not present for the 
second-order equation ( |l. 2| ) . The first is that, when n > 2, it is no longer possible to diagonalize the 
nx n matrices in the associated first-order systems in the same explicit way as when n = 2. We overcome 
this difficulty by adopting an approximate diagonalization process and this approach is discussed in 
sections 2 and 3. Second, this approximate procedure necessarily involves differential systems whose 
coefficient matrices contain a greater number of terms than when n = 2, and we have to develop a 
new algorithm for generating and collating these terms. These matters are introduced in sections 4-5. 
Then in section 6, we develop a computational algorithm to perform the repeated diagonalization. The 
estimation of error terms leading to em presents fewer difficulties and it is also dealt with in section 6. 
The application to the spectral matrix (mij(A)) is naturally less simple than when n — 2 because the 
underlying spectral theory involves several i^(0, oo) solutions of the differential equation rather than just 
one such solution. We cover this application in sections 7 and 8 with special reference to the example 
q{x) — —x" (a > 0). In section 9 we are able to give independent confirmation of our results when a — 1 
in terms of the higher-order Airy equation. Finally, in section 10 we indicate possible extensions of our 
work and, in both sections 8 and 10, we comment on the effectiveness of our methods as compared with 
the recent alternative approach of Bennewitz et al. 0. 



2 Diagonalization in differential systems 

In this section we introduce the theoretical basis for estimating and improving error terms in the solution 
of differential systems. These error terms give rise to cm in ( p. ID and in the corresponding formulae for 
( |1. 3| ) (see ( |3. 14 ) in the next section) . The main components of the discussion are covered in the three 
subsections (a)-(c) below, and we frame the discussion in terms of the system 

Z'{x) = xf^{D + 0{x-')}Z{x), (2. 1) 

where 13 is a constant diagonal n x n matrix, 

D = dg(di, ...,d„) 

with distinct dk and the O— term refers to a; — > oo. Also, /? > and 7 > 0. If in addition /3 — 7 < —1, (^ 
|l]) has the standard Levinson form section 1.3 ] 

Z' = {A + R)Z (2. 2) 



in which A is diagonal and R is L{a, 00), and the Levinson asymptotic theorem |11, Theorem 1.3.1] 
states that there are solutions Z^ {1 < k < n) of ( p. 1| ) such that 

Zkix) = (efc + rik) exp( / t'^dfedt) (2. 3) 
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where ek is the unit coordinate vector in the fc— direction and rjk ^ as x —^ oo. The size of rjk is related 
to the size of x^^^ as x — > oo in a manner to be made precise in section 3. 



(a) Exact Diagonalization We first consider the effect of expressing the coefficient matrix in ( |2. l| ) in 
its diagonal form and, as a start, we suppose that the 0-term is simpley x~^C , where C is a constant. 
Thus we write 

T-^{D + x-''C)T ^ Di, (2. 4) 

where Di is the diagonal matrix formed by the eigenvalues of D + x^^'C and T is the approximate identity 
matrix formed by the eigenvectors. Then 



Di D + Oix~'^), T = I + 0{x'"'), T' = 0{3 



-7-1 



as a; ^ OO. The transformation 



therefore takes ( |2. l| ) into 



TW 



W = xf^iDi - x-Pt-^T')W = xl^iDi + 0{x-P~'i'^^)]W. 



(2. 5) 



(2. 6) 



Hence the effect of the t ransf ormation ( 2. 5 ) is to replace the perturbation x ^ \n ( |2. l| ) by one of smaller 



magnitude x 

Even th ough ( |2. 6|) do es not have quite as simpl e a fo rm as we specified in l[ ) , the argument which 
leads from ( |2. l| ) to ( 2. 6 ) can again be applied to ( |2. 6 ) to yield a sequence of systems 



W^^xP{D.^+R-n,)Wrr 



(2. 7) 



with approximately constant diagonal matrices D^^ and perturbations i?„j of mag nitude a;-'"('3+i)-'' ( m = 
1,2, ...). The Le Vinson formula (2.3) can be applied to any system in the sequence, and the term rjk will 
have a corresponding order of magn itude. The transformation back to the original system ( |2. l| ) via the 
sequence of equations such as ( ^. 5 ) thus gives successively more accurate asymptotic representations of 
the solutions of l| ). 

The details of this process were given in sections 2-3] for the case where ( 2. 1 ) arises from (1.2) and 



all the matrices are 2x2. These details depe nded on the explicit knowledge of all necessary eigenvalues 
and eigenvectors in equations such as (|2. 4| ) . When n > 2, however, s uch e xplicit knowledge is not 
available in ( 2. 4 ). In this paper, we overcome this difficulty by modifying ( 2. 4 ) and considering instead 
an approximate diagonalization process which is based on the one introduced by Eastham ^ [ pT| section 
1.7]. 



(b) Approximate Diagonalization The matrix T in (|~|) has the form T = I + P, where P = 0{x- 
and the idea now is to modify (2. 5) to 

Z={I + P)W (2. 



with a different, but explicit P. It follows from (l-'''.2) and (1.6.13)] that (2. 8) represents an 

approximation to (2. 5) if P is defined by 



PD - DP = x-^(C - dgC) 



(2. 9) 
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with dgP = 0. Thus the entries pij in P arc defined by 



Cij/{dj~di) (i^j). 



(2. 10) 



By (2. 



and (2.9), we obtain in place of (2. 6) 

W = x^{D + x"'{I + P)-\CP + dgC)-x-'^{I + P)-^P'}W 



'{I + P)-^{CP-PdgC)-x 



-fi 



{I + P)-^P'}W. 



Here 



D, = D 



(2. 11) 



(2. 12) 



and, by ( ^. 10| ) , the other two groups of terms arc respectively 0{x~'^'^) and 0(a;^^^^^^), of which the 
latter dominates if 7 > /3+ 1. We note that ( |2. llD contains more terms than (2. 6), but this is the price 
to be paid for having explicit terms derived from ( |2. 10 ). 



(c) Further approximation Repetition of the process ( |2. 8D - (|2. Ill), but with (|~ ll| ) as the starting 



point, involves the diagonal entries of Di in place of dj and di in ( 2. 10|) . Thus powers of x appear 
increasingly in the denominators of the P— matrices. This causes problems when we seek to develop 
computing and numerical procedures for implementing the repeated transformations, particularly becaus e 
repeated differentiation of such matrices is involved, as indicated by the appearance of P' in (|2. ll| ). 
Accordingly, we shall retain the original constant entries dj and di in the denominators as we go through 
the process. This further approximation adds new terms to the coefficient matrices {...} in the sequence 
of systems corresponding to (|2. ll| ). It is this last method that we adopt in this paper and, after these 
introductory remarks, we defer further details to section 4 where we deal with the actual system ( ^. l[ ) 
which arises from ( |l. 3| ). 

3 Asymptotic formulae for solutions 

The dominant form of the solutions of (1.3) for large x is known subject to suitable conditions on Q JTl] , 
section 2.9 ]. Our aim here is to obtain an explicit estimate for the error term in the asymptotic formulae, 
which is suitable for our purposes in sections 4-6. We begin with the system formulation 

Z' = {A + R)Z (3. 1) 

of (1.3) which is given in ||ll|, section 2.9 ]. Here 

A = gi/«dg(cc.i,. ..,.;„), (3.2) 



where 

are the n— th roots of unity, and 



LOk — exp{2(fc — l)7ri/n} 

R = -n-^A^n, (3. 3) 
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where 

Ai = n-i(Q7Q)dg(0, 1, ...,71 - 1), (3. 4) 

while ft and fl^^ have (j, k) entries tofT^ and nr^uj^ ''^ respectively. The connection between Z and ?/ 
in ( |3 iD and (|l~J) is 

r = {dg(l, Qi/", Q2/", Q("-i)/")}f]Z (3. 5) 

where Y has components y,?/', ...,?/("^^^ It follows from (|3~2|) - that can be written in the 

form 

Z' = Q^/''\D + Q'Q-^-^/"C)Z, (3. 6) 

where C is constant and 

D = dg(wi,...,w„). (3.7) 

In the case where Q and its derivatives sufficiently resemble {a > 0) and its derivatives, the system 
( [3. 6 ) essentially has the form ( |2. l[ ) with 7=1 + a/n, and the ideas in section 2 for improving the error 
term apply. At this point, we note that the diagonal entries in C are all — (n — l)(2rt)_Ml^, (2.9.25) ] 
and this, together with ( 2. 12[ ), governs our definition of Dm and in ( |3. 10| ) and ( [4. 4 ) below. 
The sequence of transformations to be defined in section 4 takes (3.6) into the form 

^'m — Q^^^^{Dm + Rm)Zmj (3. 8) 

where 

M-l 

^ = { n + P'n)}ZM = {! + P)Zm, (3. 9) 



Dm = D-{n- l)i2n)-^Q'Q-^-'/"I + Am (3. 10) 

and Rm G L{X^ 00). Also P{x) and Am{x) are both o(l) as a; — > 00. We can now state and prove a lemma 
which gives the asymptotic form of the solutions of (1.3), incorporating an explicit estimate of the error 
term derived from a knowledge of Rm- In the lemma we write P = (pij) and Am = dg((5iAf , SnAi)- 

Lemma 3.1 Let M{> 2) be an integer and let Dm, Rm o-nd P he as in ( \3. - ^3. IL ). In some 
interval [X,oo), let 

RcUlu, ~LUk + S,M - 4m)Q'/"} (1 < J, fc < n) (3. 11) 

have constant sign (either > or <0), with 

poo 

I / Re{(c^, - ujk + 5jM - 4M)Q'/"}rft h oo (3. 12) 

for j ^ k. Also, let 



X 



)^/"{t)RM{t)dt \< oo. (3. 13) 

Then, for 1 < k <n and I < r < n, (1.3) has solutions such that 



X 



(r-l) 

Vk 



Q-("+'-'''^/'"K-' + E^r'(^J+Pjfe + EP^'^')}exp(/ Qi/"(t)K + 4M(t)}di) (3. 14) 
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where, in [X, oo), the rfj (1 < j < n) satisfy the estimate 

r}j{x) |< / I Qi/"(i) III RM{t) II dt/{l - n / | Q^'^{t) \\\ RM{t) \\ dt) 
with II Rm 11= max | r^jM \ {I < j,k < n) in terms of the entries in Rm ■ 



(3. 15) 



Proof. By ( |3. and (3. 13) , we can apply the Levinson Theorem |T^, Theorem 1.3.1 ] to the system 
( [3. 8 ), and then transform back to Y via ( [3. 9 ) and ( 3. 5| ) to obtain vectors Y with the form 



Y = dg(i,gi/",...,g("-i)/")r!(/ + p)(efc + u) 

X exp( r[Qi/"K + 4m} - {n - l){2n)-\Q' /Q)]{t)dt) , 



(3. 16) 



X 



where u = o(l) and we have used ( |3. 10| ). Let u have components rjj {I < j < n). Then (3. 14) follows 
from ( |3. 16| ) on taking the r— th component on each side. 

It remains to establish ( 3. 15| ) and, to do this, we proceed as in ^ (1.4.23) and (1.4.13)]. It follows 
from ( |3. ll[ ) and (3. 12) that, in the notation of ||ll|, (1.4.13) ], the non-zero entries of $i(a;)$~^(f) and 
^2ix)^~^{t) are all < 1 in modulus. Hence O, (1.4.13)] (with a = X here) gives 



Thus, a H — max | r]j \ {X < x < oo, I < j < n), we obtain 



nOO / " \ 

Jx V ^ ^ ' Q^^'^rjiMVi I 1 dt. 



H < 



-,1/r. 



R 



M 



-Hn I Q 



l/r. 



Rm II dt. 



X 



Hence 



H< I \Q 
Ix 



l/n 



II i?M ||dt/l 1-n I 



Rm \\dt\, 



and ( |3. 15| ) follows. 



4 The sequence of transformations 

In this section we define the transformations 

Zm = {! + P™)^™+i (1 < m < M - 1) 



(4. 1) 



which lead from ( 3. (: ) to ( 3. 8 ) in the manner introduced in section 2(c). A typical system in the process 
is 

=Ql/"(An+i?m)^m (4.2) 
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with (3. 6) being the case m — 1. In ( 3. (\ ), we emphasise the role of the diagonal terms in C by taking 
these terms over to D as in ( 3. 1C| ). Thus we define 



D,=D+^{n- l)pl, i?i = Q'Q-i-i/"(C - dgC) 



(4. 3) 



with D as in and p = (Q^^/")' -n-^Q'Q-^-^/". Then, again as in ( [3. lOD , we write 

A„ = Di + Ara {m > 2). (4. 4) 

Now, as indicated by ( |2. llD , Rm will contain terms of different orders of magnitude as x —> oo, and 
we wish to identify these terms according to their size. Hence we write 



where 
and 



Vm + Em — Vim + V2,n + ••• + Vnm + E„ 



Vkm = 0{Vjm) {X ^ (X),k> j) 



E„ 



o(V^m) (X -> OO). 



(4. 5) 

(4. 6) 
(4. 7) 



Here E^ represents terms which are already of the accuracy that we require in (3. 15), while the Vjm 
represent terms which are not of that accuracy and which have to be replaced by smaller-order terms as 
we go through the transformation process in the manner discussed in section 2. Also, as in the case of 
( [f. 3| ) , we arrange that 

dgl^i^ = 0. (4. 8) 

To discuss a typical step in the process leading to (3. 8), we substitute ( [4. l| ) into ( |4. 2]) to obt am 

{I + P,n)Z'm + l + PmZm+l 

(4. 9) 



Now we define Pm by 



PmD - DPm = Vir 



(4. 10) 



with dgP,„ = 0, this definition being consistent in the diagonal entries because of ( 4. 8 ). Thus the entries 
Pijm in Pm are defined by 

Pijm = Vijim/[ujj -UJi), (4. 11) 

again as indicated by ( |2. 9D , ( ^. 1C| ) and section 2(c). It follows from ( [l. 4| ) and ( [4. lOD that, in ( |4. Qj ), 
we have 

EmPm PraEm ^" ^Irn ^iriPm Em^rn Efn^ (4. 12) 



say. Hence, so far, (4. 9) gives 
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7' 



+ {Rrn - Vl„,){I + Pm))}Zra+l. 



We wish to show that this can be expressed as 



ZL 



m + l 



Q^^" (Dm+l + Rrn+l)Zm+l, 



where Rm+i has a similar form to (4. 5): 

Rm+l = Kn+1 + Em+1 = ^l,m+l 



V. 



/1, 771+1 



'm+l ' 



(4. 13) 
(4. 14) 
(4. 15) 



but with a different /i, and where the Vj^m+i can be obtained constructively from the Vjm- Also, as in 



( [1. 8 ), we shall arrange that 



dgVi,„+i = 0. 



(4. 16) 
(4. 17) 



-1 di'+It/. 

V 77 



We establish ( 4. 15| ) by writing 

{I + P„,)-' = I - P„. 
Here v is chosen so that the product 

{I + Pm) Pm 'jm, 

which occurs in ( |4. 13| ), has a sufficiently small order of magnitude to be included with and form 
part of Em+i- Thus will differ for different j, and similarly for the other terms in ( [1. 13[ ). We now 
group t ogethe r terms of the same order of magnitude and denote the dominant term by Sm+i- We then 
obt ain (4. 15 ) ( with Sm+i in place of Vi^m+i), where E^+i has the same order of magnitude as Em and, 
by ( 4. 11 ), S„i+ i and the Vj^m+i are k nown explicitly in terms of the Vjm- We complete the derivation 
of (|4. 14| ) and ( [4. 15| ), in which (|4. 16| ) holds, by defining 



and 



Dm+l — Dm + dgS'm+1 



Vl,„i+1 — Sm+1 — dgS'm+i. 



(4. 18) 



In the next two sections, we deal with the question of developing an algorithm, to be implemented in the 
symbolic algebra system Mathematica , for determining the Vjm (jn = 1,2,...) and estimating the Em- 

5 The basis of the algorithm 



The details of the procedure based on ( 4. 5 )- (t4. 18 ) depend on the nature of Q. In this section we assume 
that 

Q{x) ^ (const. )x" {x oo), (5. 1) 

where a > 0, and that differentiation can be carried out in the sense that 



)('■) 



(x) = (const.)a;"-'^{l + o(l)}. 



(5. 2) 
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At the end of the paper, we indicate the modifications which are made to accommodate other types of 



We start the detailed examination of (4. 5)-(4. 7) with the cases m = 1 and 2. When to = 1, (4. 2) 
is (B. 6) and, as in ( 4. 3| ), we have 



Di=D+^{n- l)pl, i?i = Q'g-i-i/"(C - dgC) = 0(.T-'^), 



by (|5~1) and where 

a — 1 + a/n. 

So far, we have only fi — 1 and iJi = in ( |4. 5| ). Then ( |4. 13| ) is simply 

z^ = gi/"{i^i + (/ + Pi)-i52}Z2, 

where 

S2 = -Q-^^''Pi + RiPi. 



By (|4. llD and (|~l|)-(|rj), we have 



(5. 3) 
(5. 4) 

(5. 5) 
(5. 6) 

(5. 7) 



For the inverse matrix in ( 5. 5| ), we use ( |4. 17| ) and, at this point, we specify the accuracy to be 
represented by the Em in (4. 5). We require 



(5. 8) 



for all TO and a fixed integer > 0. Thus the process leading to ( 3. 8 ) ends when ( 4. 5 ) reduces to 

Rm = Em. (5. 9) 

Then by ( |4. 171) , we have 

(/ + Pi)-'S2 = (^E(-l)-'"'^l ^2 + + Pl)-'Pl'S2. 



Thus (|5. 5| ) is the case to = 2 of ( [4. 21) , where ( |4. 5]) holds with n = N, Vj2 = (-l)^^^Fi'"^S'2 
(2 < j < iV) and 



Also, as in (4. 18), 



We note that, by ( |5. 7] ) 



and, by (4. 4) 



ii;2 = (-1)^(7 + Pi)-i<52 



D2 = Di+ AgS2, V12 = 82- dgS'2. 

v,2 = o(x-(-''+i)'^), i;2 = o(x-(^+2)'^) 

A2 = dg^2 = Oix-^''). 



(5. 10) 



Then, since Am (to > 3) is obtained from { i. 18[ ) by adding successively smaller-order terms to A2, we 
have 

A™ = 0(x-2-) (5. 11) 

for all TO. 

We move o n to g eneral to in ( |4. 5 ). An easy induction argument starting from ( 5. 3| ) and ( |5. 10| ), 
and based on ( 1. 11 ) and (|4. 17| ), shows that 



0{x 



-(m+j-l)a 



Then, by (|5. l| ) and (^. ll| ), it follows from ([4. ll| ) and ( [4. 12|) that 

Pm = 0(x-""), Q"'/"P;„ = 0(x-("+i)"), Tm = 0(x-("+2)°). 
It now follows from ( |4. Sj) , ( |5. 121 ) and (|5. 13| ) that, in ( fl. 14]) and ( [4. 15| ), 



(5. 12) 



(5. 13) 



(5. 14) 



and then (5. 14) is used in ( |4. ll| ) to define Pm+i- A further consequence of ( [4. 5D and ( ^. 12| ) is that 
(|~9) is achieved when 

M = iV + 2. (5. 15) 

It is not difficult to check that the order relations ( 5. 12| ) and ( 5. 13| ) are in fact exact, and we can 
describe ( |5. 12| ) by saying that Vjm has exact order m + j — I, and similarly for ( ^. 13 ). The grouping 
together of terms of the same order is the basis of the algorithm which we go on to describe in more 
detail now. 



6 A symbolic algorithm for the computation of solutions 

We can now use the ideas i n sec tion 5 to discuss the development of a co m putat ional algorithm for 
computing the solutions of (|l. 3| ) over an interval [X, 00), subject to ( |5. l[) ( ^. 2| ). The algorithm is 
implemented in the symbolic algebra system Mathematica and we use the notation in sections 3-5 to 
denote the symbolic objects that we need. The algorithm computes estimates for n linearly independent 
solutions of ( |l. 3| ), and for the derivatives, and it is structured in three distinct stages. In the first stage, 
we shall not make any assumptions about the order of the differential equation while, for reasons of 
clarity, in the second and third stages our discussion is focussed on the case n = 4 and more particularly 
on the equation 

y'-'^^ -x"y = Xy {X < x < 00), (6. 1) 



where X > and a > 0, this equation being covered by ( |3. l| ) and (5. 2). 



We recall that the discussion in sections 4-5 of this paper shows how the system 

Z;„ = Oi/"(i?„, + Vra + E„,)Z,„ (6. 2) 

is transformed into 

^m+l ~ Q^^"(-Dm + 1 + Vrn + 1 + -E'm+l)^m+l (6. 3) 
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by the mapping 

Zm = [I + Pm)Z^m+l. (6.4) 

Here P„i, D„i and Vm are defined in ( [4. 11| ), ( [4. 18| ),( [4. E| ) and (4. 8), where Vjm has exact order m+j — 1 
as described at the end of section 5. We now restate some of the features of this discussion in a suitable 

form for implementation in our algorithm. 

First, the number M which appears in (B. 8) and (5. £) is determined by the accuracy of our working 
which, by (|7J) and (|5. ISj ), can now be stated as 



E„ 



Oix-'''") (2 < TO < M). 



(6. 5) 



Then, by (0|), we have /i = M__rn in ( [4. 5| ). Next, examining the use of (g[T|) in (ITTsI ), we denote 
by U any one of the terms in (4. 13) on which the inverse acts, excepting E„i- Thus 



U e {-Q-^^''P^,,Tm, VlmPm, Vjrn, V.mPrn (2 < J < M ~ m)} . 



(6. 6) 



For each [/, the integer v in ( |4. ITf ) is chosen so that P'^^ U has sufficiently small order of magnitude 
to satisfy (6. 5). Since the order of U is known from (5. 12) and (5. 13), the computation of v and the 
grouping of terms {—lyP^U {Q < r < v) oi the same order can be performed for each U and r. 
By (4. 5), the error term Em appears in ( 4. 13[ ) in the form 

AmEm{I + Pro), 

where the symbol stands for (/ + Pm)~^ ■ This is taken as an initial definition of E^+i which is then 
updated as ( [4. 17| ) is used for each U, to yield the final E^+i in ( 4. 15\ ). 

This discussion leads to the following algorithm, based on (4. 11) and ( |4. 17D , for computing the Vjm, 
and Sm — and hence the Dm and Vm — in sections 4 and 5. 

Algorithm 6.1 ( ^ ) First input M to fix both the number of iterations required for ( jj'. ^ 
and the order (6. 5) of Em- 

( 2 ) Start with Di and Vi, and put Ei = 0. 

(3)Form = ltoM-l, Em+i - AmEm{I + Pm)- 

( 4 ) Por each U in (6. 6), determine v, the number needed in 1'^ . 

( 5 ) Update error term, i.e. Em+i — Em+i + i^^)'^ AmP^^^U . 

( 6 ) For r — Q to V , determine the exact order p — mr + (order of U) of PmU . 

(7) Updat e V,-m .m+i ^ Vp-m.m+i + {-If P^U {p-m>2). 
Then, as in ( 5. 14 ), 

Sm+1 = —Q '^^^P'm + ^2m 



and 



Vi^m+l — Sm+1 — dgS'm+l. 
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At any stage in the algorithm, ^m+i depends on the terms Dj,Pj,P^ and Vj (1 < j < m). However, 
simpHfications can be made by retaining earHer terms Sj where possible. This reduces the number of 
terms in the computation of Sm+i with consequent economy in the algorithm. We illustrate this point by 
giving the results for S3 and S'4 which are first obtained in terms of Vi, V2, Pi etc, and then the previously 
computed definition of S2 is used dynamically to simplify the expressions for S3 and S'4. Thus we have 

52 = + (as in (5.6)) 

P2 

53 = —-—Y + Ti~PiS2 

Si = -^+T2-PlTi + V2P2+P?S2. 



Q 

We also give the error term as an example 



= AU3A2Ai{Pi^Ti-Pi^S2}iI + P2)iI + P3)iI + Pi) 

+ { (-P1T1 + Pi^S2 - P1S2 + Ti) P2 

-P2 {S2 - PiTi + Pi^S2 + T2 + V2P2) } (/ + P3W + Pi) 
+ AiA3 {T3 - P3S2 + V3P3 + (-P1T1 + Pi^S2 + V2P2 + T2) P3 } (/ + Pi) 

+ AA-^+Ti + ViPA. 



The expression 

Se = -Q-^P;,+Ti+{-PiTi+P^S2)P2 

- P2 (-P1T1 + V2P2 + P?S2 + T2) - PfTi + P^S2 + V3P3 

which we use later in our examples, has an associated error term Eq with over 250 components. The 
numbers of components in Sg and E^ are about 60 and 700 respectively. This complexity indicates that 
the computation is intractable by hand, and that a symbolic algebra system is the only way to generate 
the recurrences. The algorithm is however quite cheap to compute, and we give some sample times in 
Table % 



m 


duration in seconds 


4 


1.6167 


6 


7.5000 


7 


22.4167 


8 


61.9833 



Table 1: Time in seconds for computing expressions Sm and E, 
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At this stage, the algorithm involves nothing more than a series of recurrence formulae in terms of 
objects which satisfy non-commutative multiplication and certain order relations. When we wish to go 
further and apply the algorithm to a specific equation ( |l. 3| ) , we need to express the recurrence formulae 
in terms of n x n matrices and, in the next stage of the algorithm, we focus on n = 4. In this case (3. 7) 
becomes 

Z? = dg(l,i,-l,-z) 



and, in ( 4. 3 ) 



Algorithm 6.2 Starting with initial matrices 



(6. 7) 



Also, it is easily verified from ( 3. 3 ) and ( 4. 3 ) that Vi (= Ri) has the form to be stated now. 



Vi 



Di = 


D + 


2^^' 







1 + 


i 1 


1 - i 


1 - i 





1 + i 


1 


1 


1 - 


i 


1 + i 


1 + i 


1 


l-i 






and El — 0, the expressions S2 

then used in turn to generate the 4x4 matrices D2 



Sm-1 generated in Algorithm 6.1 are evaluated in order. These are 

Pm-i, Pi,- 



,Dm-1 ,V2, ...jVm-I, P2,. 



p' 



We note that the process stops at M — 1 because ( ^. 9 ) implies that Dm = Dm-i and Vm = 0. In 



Algorithm 3.2, we assume only that p and Q are sufficiently smooth functions related by ( ^. 7] ) and 

Q" = -AQip' + 2Qp^Qi. (6. 8) 

This equation is used to express higher order derivatives of Q in terms of those of p. 

There are several severe computational problems contained in this seemingly simple process of matrix 
multiplication and substitution. Although all the operations needed are to be found in the Mathematica 
system, the use of the Substitute command in Mathematica requires an inordinate amount of computa- 
tional effort, and therefore it is necessary to consider what substitutions and simplifications are best used 

have been found most useful in sim- 



at each stage of the algorithm. The relationships in (6. 7) and (|f 
plifying the expressions produced by this part of the algorithm. The main problem however, manifested 
by long computing times, lies in the exponential growth of the number of terms in the elements of the 
matrices Sm- For example, each of the elements of Sf, consists of a sum of over one thousand terms, each 
of which is a product of several distinct objects. A consequence of this complexity is that, as m increases, 
the amount of computational effort required to perform the calculations increases exponentially. The 
time in cpu seconds needed to compute Si and S^ on a SPARC 10 workstation, is given in Table ^. The 
- indicates that is has not been possible to complete the computation within a reasonable time. This is 
because of the limitations in the size of memory and cpu speed that we have available. 

The final stage of the algorithm deals with the remaining matter of obtaining an explicit upper bound 
f or the norm || || of E^n which adds precision to the order estimate ( |6. 5| ) and which can be used in 
( [3. 15 ) when m = M . We use the sup norm, for which the triangle and Cauchy inequalities are 



A + B\\<\\ A 



B 



AB ||< 4 II A 



B 
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TO 


duration in seconds 


4 


74.6167 


6 


904.24 


7 





Table 2: Time in seconds for computing matrices S, 



in our case of 4 x 4 matrices. For the inverse matrix A„i = (I + Pm) ^, we have from the geometric series 

II ||<i + f;i|p;;ii<i + f;4'-i i|p„, ir 

r— 1 r— 1 

and hence 



Am \\< 1+ II Pm II /(I -4 II P„ 



(6. 9) 



provided that || P^ \\< 1/4. 

There is one other point to mention before discussing this final stage of the algorithm, and this is 
prompted by the fact that the reciprocal of Q occurs throughout the process in Algorithms 6.1 and |6.2| . 
In order to estimate this reciprocal, we require a lower bound for Q in the relevant a;— interval [X,oo). 
Thus, in addition to ( ^. l|) and ( 5. 2 ), we require that 



I g(x) |> fc^^x" {X<x<oo) 



(6. 10) 



for some constant fc > 1. In (|6. l| ) itself, we have Q{x) = X + x", and we shall consider A in the half 
plane 

ReA > -1. (6. 11) 



Thus I Q{x) \>\ ReA + a;" |> (1 - X-")x°' provided that X > I. Thus ( |6. lOj ) holds with 

k = {i-x-'^yK 



Algorithm 6.3 Compute the sup. norm of each matrix in Em, using ^6. !\ ) for the inverse matrices. 
Next apply the triangle and Cauchy inequalities to obtain an upper bound for the sup. norm of Em itself. 

The terms that we encounter in computing the matrix norms involve p and its derivatives. In order 
to obtain an upper bound for the k^^ derivative of p{x), we first use 

Q{x) = A + x" 



and (3. 8) to compute a symbolic expression for p^^\ We then use the inequality ( |6. lOj ) to compute 



bounds for | p^'^^x) \. In obtaining this estimate every element of every matrix which occurs in the error 
term must be examined. The triangle and Cauchy inequalities are applied to each element and the 
bounds for I Q I and | p'*"^ | are substituted with fc > (1 — A^")^^ and the specific values of X and a 
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under consideration. This enables the sup. norm of each matrix to be computed at the point X and a 
precise upper bound e{X) for the error in the solutions determined. Again the amount of computational 
effort needed is exponentially increasing with the requested number of iterations. In table ^ we show, 
where possible, the amount of cpu time needed to compute these norm estimates. 



m 


durations in seconds 


Sm I Rva 1 Pm i UOrmS 


T,„ norms 


P/„ norms 


4 


167.833 


148.25 


241.33 


6 


1159.5 


1446.2 


2216.25 


7 


3297.95 







Table 3: Time in seconds for computing norms where k = ^ and a = 1 



With the completion of the three stages of our algorithm, we have all the information needed to 
compute the product 

m— 1 

which appears in ( 3. 9 ) and which c ontrib utes the terms pij to ( |3. 14|) . Also, by ( 5. 9 ) 
information needed for the estimate (3. 15) of rjj which represents the error term in (3. 14). Thus our 
algorithm computes the solutions ( |3. 14 ) of ( 1. 3 ) to a specifie d acc uracy based on ( 3. 15 ). We have of 



course chosen to orientate our discussion towards the example ( |6. l| ) 



(6. 12) 



we have the 



The computation associated with the product ( 6. 12| ) can be implemented symbolically only for 



M < 6 since, for larger M, the number of terms that need to be manipulated becomes too large for the 
computer store that we have available, and the computation is performed numerically with an accuracy 
of 30 decimal digits. In table ^ we give the times needed to compute ( [3. 14 ) using symbolic and numerical 
methods. The shorter time needed when a = 1 is a consequence of derivative terms, which are not zero 
for non- integer a, becoming zero when a ^ 1. 



M 


symbolic calculation 


numeric calculation 


a = ^ 


a = 1 




a = ^ 


a = 1 


a = ^ 


4 


120.65 


17.5167 


149.983 


9.55 


6.46667 


13.8833 


6 




182.583 




34.1333 


11.2 


40.5333 


7 




773.8 




71.556 


19.35 


84.6667 



Table 4: Time in seconds for determining (3. 14) with X = 10, and A = 1 
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7 Generalised Titchmarsh-Weyl theory 



The motivation for the asymptotic analysis and associated algorithm in sections 2-6 is provided by the 
spectral theory of the equation 

i-l^y'-^''^ +q{x)y = Xy {0 < x < oo) (7.1) 

together with boundary conditions at a:; = 0. Here q{x) is real-valued and locally integrable in [0, oo) and 
A is a complex spectral parameter. We take the boundary conditions to be the Dirichlet conditions 

y('=-i)(0) = (1 < fc < I/) (7. 2) 

but we shall comment briefly on other choices of boundary conditions later in this section. Our purpose 
here is to state what we need from the spectral theory of ( |7. l[ ), and then to introduce the contribution 
that our methods make. 

The fundamental result of Titchmarsh and Weyl in ||2^, Chapter 2] and |Q concerns the case v = 1 
of ( [7. l| ) and ( [7. 2| ), and they showed that, when ImA ^ 0, ( [7. l| ) has a non-trivial solution ■!/)(a;, A) 
which is i^(0, oo). The importance of this solution lies in writing in the form 

^'(x, A) = e{x, A) + m{X)(l){x, A), (7. 3) 

where 6 and ^ are the solutions of (|7. l| ) ( with = 1) which satisfy the initial conditions 

61(0, A) = 1, 6i'(0,A)=0; 0(0, A) =0, 0'(O,A) = 1. (7.4) 

The function m{\) thus defined is the basis for the spectral theory of ( |7. l| ) and ( [7. 2D as developed in 
pSf Chapter 2] and |^ for the case v = 1. Depending on the nature of q{x), either m(A) is uniquely 
determined for every non-real A ( the so-called limit-point case ) or m(A) involves an additional parameter 
( the limit-circle case). We refer to the paper by Fulton |jl^ for a more recent critical discussion of the 
li mit-circle case. A number of explicit examples of m(A) are given in |2^, Chapter 4] for equations ( 

I l|) which can be solved in terms of the special functions. 

The generalisation of the Titchmarsh-Weyl theory to higher-order equations, of which ( [7. l| ) is an 
instance, was made by Everitt , |lj] . Here we indicate this generalisation in the yet wider context of 
the Hamiltonian system 

JY'{x) = {\A{x) + B{x)}Y(x) (0 < X < oo), (7. 5) 

the spectral theory of which was initiated by Atkinson p]. C hapt er 9] and further developed by Hinton 
and Shaw in a series of papers which include ||l^, ||l^, |[l9[|. In (7. 5), F is a 2i^-component vector, A 
and B are Hermitian matrices with ^ > 0, and 



J 



-h 
h 



in terms of the v x v identity matrix ly . We note that ( [7. iD is the special case of ( |7. 5| ) in which the 



entries a.y and hij of A and B are 



Oil = 1, a.y = otherwise, 
bii = -q, b2uau = 1, 
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bij — 1 {j = i + i'— I, 2<i<i' and i ^ j + v — 1, 2 < j < v), 
and bij — otherwise. The components yi of Y are then given by 

We now introduce the Dirichlet boundary condition 

(/. o)r(o) = 



(7. 6) 



(7. 7) 



for ( 7. 5 ), which corresponds to ( 7. 2 ), and we also introduce the 2iy x v solution matrices Q{x,X) and 
$(x, A) of ( 7. 5 ) which satisfy the initial conditions 



8(0, A) = 







, <I'(0,A) = 







(7. 



These conditions correspond to ( [7. 4|) . In the references just cited, it is shown that, when ImA ^ 0, (0 
I 5[ ) has a X ;y solution matrix, with rank such that 

(7. 9) 



(7. 10) 



'^*{x,X)A{x)^{x,X)dx < oo. 

Further, corresponding to ( [7. 3| ), ^' can be written as 

xif{x, A) = Q{x, A) + 4>(a;, A)M(A), 



so defining the v x v matrix Af(A). This matrix again forms the basis for the spectral theory of (7. 5) 

and ( [7. 7| ), as in the case of to(A) in ( 7. 3 ). 

In analogy with the Titchmarsh-Weyl theory, we say that (7. 5) is in the limit-point case if, for some 
A with ImA ^0, (7. 5) has no more that linearly independent solutions such that 

Y*{x,X)A{x)Y{x,X)dx < oo. 

This limit-point classification is independent of A and, by (7. S) and ( [7. ICj ), it has the implication that 
M{ X) is uniquely determined for all non-real A. Further, if ^'i(a;. A) is any other 2^ x v solution matrix 
of (7. 5) which has rank i/ and also satisfies (7. 9), then 

«'i(a;,A) = *(x,A)C(A) (7. 11) 

for some non-singular v x i/ matrix C. It now follows from ([7~8|) , ( |Oo| ) and ( fTTll ) that 



«'i(0,A) 



Hence, partitioning 



*i(0,A) 



C(A) 
Af(A)C(A)) 



C(A) 
77(A) 



(7. 12) 
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we obtain 

C(A) = C{X), r,{X) = M(A)C(A). 
Eliminating C(A), we obtain a formulae for M(A) in terms of the initial values of "ifi a,t x — 0: 

M(A)=,7(A)r^(A). (7.13) 

Finally here, we note that M(A) is analytic for ImA 7^ and M(A) coincides with its own transpose, that 

is, its entries satisfy niij(X) — mji{X) [ |l8| . 

For use in the next section, we mention here the form which ( [7. 13|) takes when (|7~|) arises from the 
fourth-order example of ( [7. l| ) with ly = 2. In this situ ation , w e have two li nearly independent L^(0,oo) 

solutions ■ipi{x,X) and i/'2(a;, A) of ( [7. l|) . Then, by (|7. 6| ), ( |7. 12| ) and (7^ 13), the spectral matrix 

M{X) = (my (A)) is given by 



,(A)) 



-<'(0,A) 



-V^^"(0,A) 



V^1'(0,A) ^^'(0,A) 



^1(0, A) 
^i(0,A) 



^2(0, A) 
^2(0, A) 



(7. 14) 



We have mentioned that examples of m(A) exist when i/ = 1 in ( |7. l| ), arising from the special 
functions. In contrast, no such examples are known when v > 2 except for the relatively trivial Fourier 
case when q{x) = 0. In order therefore to illustrate and test the general spectral theory of ( |7. l| ), there is 
a need for a computational approach to estimate M(A) numerically. Thus the motivation for our work is 
the development of an effective computational approach and, in the next section, we give our numerical 
results for the example where v = 2 and q[x) = — x". Our methods are not confined to this choice of 
g(x), but this choice does have a certain significance which we discuss in section 10. 

We have concentrated on the Dirichlet boundary conditions ( |7. 2[ ) and (7. 7). If these are changed 



to other conditions which involve linear combinations of derivatives, then (7. 8) and consequently M(A) 



are also changed. However, in the limit-point case emphasised here, the new M(A) is related to the 
Dirichlet M£)(A) by a functional equation, and therefore in principle the properties of the new M(A) 
can be obtained from those of My){X). We refer to p5|, p. 69] and P, Theorem 3.4] for the details of 



the functional equation. Here, as an example, we note only that, if (T 
condition (0/^)^(0) = 0, then §, Corollary 3.1] 

Mn{X) = -M^^{X). 



7) is replaced by the Neumann 



8 Computation of the spectral matrix 

We can now draw together the ideas in the previous sections to give an effective procedure for the 
numerical computation of the spectral matrix (mij(A)) associated with the equation 



V 



(4) 



x°'y = Xy (0 < X < 00) 



1) 



and the Dirichlet boundary conditions 



y(0) = 0, 2/'(0) = 0. 
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The range for a is < a < 4/3, so that ( |8. l| ) is in the hniit-point case as required in section 7 p] , 
Theorem 3.11.1(b)] and ||l|, section 23.4]. 



The mij{X) are given by (7. 14) and we therefore require the values at a; = of two L^(0, oo) solutions 



■01 (x, A) and ■02(3^, A) of (8. 1), together with the values of their derivatives. We take ImA > to be 



definite. Then consideration of the exponential factor in (3. 14), with Q = A + x° and n — A, shows that 
the two L^(0, oo) solutions of ( 8. 1 ) are associated with uj2 and Thus ( 3. 14) pro vides the values of 



ipi(x,\) and ip2{x,X) and their derivatives at x = X, with an error given by ( 3. 15 ). These values are 
produced by our algorithm as explained at the end of section 6, and we denote the maximum error in 
this computation by £{X). The values of the solutions and their derivatives at X are used as initial data 
for a numerical initial value solver which is then used to compute numerically the values of the solutions 



and their derivatives at 0. These values are then used in ( 7. 14 ). This computation back to a; = has 
been performed with the NAG library code D02NMF in the following way. We introduce the standard 
Riccati variable ^{x) = a{x)T^^ (x), where 



Then, as in ^ satisfies a first-order non-linear differential equation. The above code is applied to ^, 
and it computes £, back from X to 0. Then, by ( [7. 14] ), Af (A) = C"^(0). 

The value of X is at our disposal, and the actual choice is governed by two factors. A smaller choice 
has the advantage that the integration back from X to zero is performed over a smaller range with a better 
retention of the accuracy represented by e(Ar). On the other hand, a larger choice has the advantage 
that e(^) is smaller. The first factor appears to be more significant and we have chosen X = 10. Later 
we shall comment on X=20 and on the very large values of X that arise in the alternative method of 
Bennewitz et al . The value of M — the number of requested transformations in sections 3 to 6 — is also 
at our disposal, and we have chosen M = 6 to obtain the good accuracy which we present later in Tables 
5 We turn now to discuss the results of numerical experiments performed with our algorithm. In Tables 
5 -0 we give the values of toii(A), m22(A) and mi2(A) = TO2i(A) for three values of a and a selection of 
values of A. The values of a are, first, a = 1 because in this case we can obtain independent confirmation 
of our results from the theory of the higher-order Airy equation, which is the subject of section 9. Then 



a = 4/3 is chosen because it is the maximum value of a for which (|8. l| ) remains in the limit-point case. 
Finally a = 1/2 is chosen in order to give a spread of values of a in the range < a < 4/3. 

The values of A, all with ImA > 0, are mainly taken to be near the origin because we wish to compare 
the power of our method with that of the code developed previously by Bennewitz et al. ||^. Our 
method is considerably more effective, and the improvement is most conspicuous for smaller values of 
A I, particularly when a = 4/3, as can be seen from Table ^ and the choice of X in the other tables. We 
recall that a value of X as small as possible is generally the most desirable because of the integration of 
l\ ) over [0,X] starting from initial values at X. 

8.1 Example a = 1 

We have carried out the computation using Af = 6 in the algorithm in section 6. With X — 10 and 
X = 20, we obtain 

e(10) = 2.036696 x 10~^ e(20) = 1.038152 x 10"*. 
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\ 


Asymptotic solution 


Bennewitz 




numerical integration 
over [0, 10] 


over [0, 20] 


accuracy 


integration 
distance X 


i 


-0.1844321489 + 1.8220405579 i 


10 sf 


10 sf 


11755.07 




-0.6613801122 + 1.1030735970 i 


10 sf 


10 sf 






-1.4291142225 + 0.8401392698 i 


10 sf 


10 sf 




0.5 + i 


0.1935073882 + 2.0560526848 i 


8 sf 


Failed 


74606.02 




-0.5871452689 + 1.2561461926 i 


8 sf 








-1.4407042265 + 0.9327853918 i 


8 sf 






10 + 10 i 


2.0868480206 + 10.2971620560 i 


10 


10 sf 


59.15 




-1.4532939196 + 3.5441000462 i 


10 


10 sf 






-2.3046081066 + 1.5476453304 i 


10 


10 sf 




1 + 0.001 i 


1.0726006031 + 1.6053582430 i 


7 sf 




Failed 




-0.2162801623 + 1.3155398369 i 


7 sf 








-1.2967267036 + 1.0783087015 i 


7 sf 






0.01 i 


0.3271814883 + 1.0330723524 i 


7 sf 




Failed 




-1.2149009705 + 0.8802024126 i 


10 sf 








-0.3127440214 + 0.9515559077 i 


7 sf 







Table 5: a = l 



The numerical integration over [0,X] is performed with an accuracy of 10~^°. Eight values of A are 
considered, as listed in Table ^. In the third column, the agreement in terms of significant figures is given 
for the choice X = 20. The slight loss of accuracy is in line with the comments which we made earlier in 
this section concerning the choice of X. The last two columns give the corresponding performance of the 
code of Bennewitz et al ||^. The large values of X and the cases of failure will be noted, as will also the 
cases of agreement. 

8.2 Example a = 4/3 

We have again carried out the computations using AI — 6, and we obtain 

e(10) = 6.57898 x 10"^ 

Table I gives the values of for eight values of A as before and the choice X — 10. In the last two 
columns of the table, the poor performance of the previous code |^ will be noted. 

8.3 Examples = 1/2 

Again with M — 6, we obtain 

e(10) = 1.452392 x 10"^ 



20 



A 


Asymptotic solution 


Bennewitz 






accuracy 


integration 
distance X 


i 


-0.2790325582 + 1.8323447704 i 




Failed 




-0.7153936028 + 1.1025729179 i 








-1.4712435007 + 0.8293465972 i 






0.5 + i 


0.0974698886 + 2.0557093620 i 




Failed 




-0.6384261250 + 1.2467917204 i 








-1.4768083096 + 0.9157238603 i 






10 + 10 i 


2.0430564880 + 10.2711343765 i 


6 sf 


290.36 




-1.4629243612 + 3.5318126678 i 


6 sf 






-2.3070857525 + 1.5415457487 i 


6 sf 




1 + 0.001 i 


0.9584534168 + 1.5867842436 i 




Failed 




-0.2741018832 + 1.2855821848 i 








-1.3294038773 + 1.0418083668 i 






0.01 i 


0.2010058761 + 1.0459299088 i 




Failed 




-1.2738007307 + 0.8492208123 i 








-0.3926285803 + 0.9405522943 i 







Table 6: a = | 
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\ 

A 




Bennewitz 






accuracy 


Integration 
distance X 


i 


0.0091688652 + 1.8079411983 i 


7 sf 


174.35 




-0.5696548223 + 1.1018460989 i 


7sf 






-1.3599216938 -1- 0.8523176312 i 


10 sf 




0.5 + i 


0.3899112046 + 2.0596635342 i 


6 sf 


193.55 




-0.4997293651 -1- 1.2700277567 1 


7 sf 






-1.3820117712 -1- 0.9563996792 1 


7 sf 




10 + 10 i 


2.2008070946 -f 10.3476095200 1 


6 sf 


33.84 




-1.4308696985 + 3.5664336681 1 


6 sf 






-2.2991588116 + 1.5578293800 1 


7 sf 




1 + 0.001 i 


1.2971330881 + 1.6384627819 1 




Failed 




-0.1199851707 + 1.3626530170 1 








-1.2458260059 + 1.1335541010 1 






0.01 i 


0.5758281350 + 1.0167049170 1 




243897.42 




-1.1166948080 -1- 0.9334350824 1 








-0.1757378578 + 0.9719497561 1 







Table 7: a = i 



and the values of the rriij are given in Table 0. This time the previous code performs better but it still 
does not match ours. It will be noted, in all three numerical examples that we discuss, that our methods 
produce results that are in agreement with the results in Further the algorithm reported on in 
fails to compute the spectral matrix for values of A close to the real line. Our algorithm has no difficulty 
computing the spectral matrix at these values. 



9 The higher-order Airy equation 

Independent confirmation of our numerical results for the case a = 1 in Example 8.1 is provided by the 
theory of the equation 



which, when n = 2, becomes the well-known Airy equation. We note that (9. 1) is itself a special case of 
the equation 
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with m rational, which has been extensively investigated by Turrittin Heading ||T^, and others ( see 
also Paris and Wood pp. 188-190]). Here z is a complex variable and, when 

z = A + a; (9. 2) 

and n = 4, ( 

becomes 

y^^\x) = iX + x)yix), (9.3) 
which is the case a = 1 of l| ). Thus we consider the equation 

g = (S..4) 

and our primary source for the nature of the solutions is the paper by Heading |l^, sections 1-5]. 

We note first that a straightforward power series substitution for y in (9. 4) gives the general solution 
in the form 

3 

y(z)=^c,zVrW, (9.5) 



r=0 



where the Cr are constants and 



^ ^+5!" +5!10!" +5!10!15!" + - 

^ , , ^ 2! . 2!7! 10 2!7!12! 

= 1 + r +6!11!^ +6!TT!16!^ 

^ . . ^ 3! 5 3!8! 10 3!8!13! 15 

= ^+7!" +7!12!" +7!T2!T7!^ 
41 4101 4191141 

■'^^ 8! 8!13! 8!13!18! ' ^ ' 



and the power series converge for all z. As mentioned by Heading |16, p. 405], the fr can if necessary be 
expressed in terms of the generalised hypergeometric function 0^3: for example 

fo{z) - 0^^3(4/5,3/5, 2/5; zVS^). 

A particular choice of the constants Cr is identified by Heading as producing solutions which have 
specific asymptotic behaviours for large \ z \. The choice (37)] is 

-i)r(l)r(|)5-" 
-|,r,-i)r(i)5-« 

-^)r(-|)j(-i)5-"/^ (9, 7) 



Co 


- r( 


Cl 


= r( 


C2 


= r( 


C3 


= r( 
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Then, as in (24), (37), (38) ], we denote the solution (9. 5) with this choice by J{z). It is shown in 
@ (39)] that 

J(z) - (const.)z-3/8exp(--z5/4) (g. g) 

5 

for large | z \ with arg z tt, and wc can use this asymptotic formula to find the two L^{0, oo) solutions 

ipi{x,X) and ^12(3:, A) of ( 9. 5 ) whe n Im A 7^ 0. 

First, with z = a; + A as in ( 9. 2 ), it is clear from ( 9. 8 ) that J{x + X) is exponentially small as a; ^ cx) 
and hence is i^(0, 00). Thus 



V'i(a:,A) = J{x + X) =^Crix + Xyfrix + X), 



(9. 9) 



r=0 



where the fr and Cr are defined by (9. 6) and (9. 7). 



Second, as observed in |l6|, (25)], y{wz) is also a solution of ( 9. 4 ) whenever y{z) is a solution and 



— 1. We therefore consider 



Mx, A) = J{{x + A)e-2W5| ^ J2 c,e-2"^'^/^(a; + Xy fAx + A). 

r=0 



(9. 10) 



For this solution, the exponential factor in (9. S) is 



exp{^«(a; + A)^} = exp{^ia;' (1 + ^Ax ^ + 0{x ^))}, 



and the modulus of this factor is 



{l + 0{x *)}exp{-(ImA)a;4}. 

Hence, when ImA > 0, '4'2(x, A) is also exponentially small as a; — > 00 and is therefore L^{0, oo). 

The initial values of V'l and -02 at x = are now known in terms of A from (3. 9) and ( |9. 10| ), and 
therefore (|7. 14| ) provides an explicit formula for il/(A) in terms of power series involving A. In using 
this formula to verify our computations in Example 8. 1, we have approximated ipi{x^ A) and "02(2^, A) by 
using the first 20 terms of the power series given in ( |9. 5| ) and ( |9. 6| ). Expressions for the derivatives 
have been calculated symbolically, and the evaluation of -01(0, A), -02(0, A) and the derivatives has been 
performed with 30 digits of accuracy. 



10 Concluding remarks 

( a ) The algorithm of Bennewitz et al. jQ. We have already referred to some numerical 
results of the algorithm Q in section 8 when q{x) = —a;". A more general comment on 
the algorithm is that it works well in cases where A/(A) is meromorphic, as for example 
when q{x) — x" (a > 0). However, the computational complexity of the method causes the 
performance to degrade and even fail for examples of ( |7. iD when argA is small and M(A) 
is not meromorphic. It is partly for this reason that we have concentrated on q{x) — — x". 
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(0 < a < 4/3), the spectrum in this case fiUing the whole real axis section 24.4] with, 
consequently, M(A) being non-meromorphic. This is the wider context within which to place 
our comments in section 8 concerning the effectiveness of our methods as compared to those 
of @. 

HELP inequalities. In the case v — 1 of (|7. l| ), the behaviour of m(A) in the neighbourhood 
of A = determines the validity of what is known as the HELP ( Hardy-Everitt-Littlewood- 
Polya) inequality. We refer to [|| and Q for surveys and an extensive bibliography concerning 
the inequality. An extension to fourth-order differentia l equ ations was given by Russell [p3[ 
[^ , but recently a more systematic development for ( |7. l| ) with general v has been given 
by Dias ^ . Again the validity of an inequality of the HELP type depends on the behaviour 
of M(A) near to A = 0. Having in this paper developed a computational algorithm which is 
effective for small | A |, we propose to investigate further the application to the validity of 
higher-order HELP inequalities. 

CoefRcients v^rith an oscillatory factor. It is pointed out in Example 2.4.1] that the 
method of repeated diagonalization for ( |l. 2| ) works not only for coefficients of the type (|]] 
J_) but also for coefficients such as 

Qix) = xyix") 

with (3 > and p{t) periodic in t and nowhere zero. It seems likely that, with suitable 
modifications to the algorithm in section 6, the methods of this paper will cover such a 
coefficient in ( |1 . 3| ) . It would appear that a higher value of M than (5. 15) is needed to 



achieve the same accuracy ( p. 8 ) as before and that the grouping of terms in the algorithm 



of section 5 and 6 depends on the value of /?. These details are another matter for further 
investigation. 

Hamiltonian systems in general. We have concentrated in this paper on the equations 
( |1. 3| ) and, for the spectral theory, ( |7. 1[ ). However, our methods are in principle appli- 
cable to other differential equations and in deed Hamiltonian syste ms w hich, after an initial 
transformation, can be written in the form ( |3. l|) with A similar to (|3. 2| ). Such a A would be 

A = p{x)dg{di, ...,dn), 

where there is a factor p{x) and the dk are distinct non-zero constants. More difficult 
however — and this is the point of this subsection — are systems ( |3. l| ) where the diagonal 
terms in A have different orders of magnitude as a; — > oo. Such a situation can occur for 
example with the equation 

y^^Hx) + {P{x)y'{x)y + Q{x)y{x) = 0, 

where the middle coefficient is dominant in the sense that Q = o{P^) as x oo |l^, section 
3.5]. Thus a further stage in the development of our asymptotic analysis and the associ- 
ated algorithm would be to cope with differential equations, and more generally Hamiltonian 
systems, where this type of A occurs. 
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( e ) Automatic differentiation The algorithm discussed in this paper has been developed us- 
ing both Mathematica and Fortran??. This has resulted in some computationally intensive 
symbolic calculations as is demonstrated by the results in Tables ^ and ||. A possible alter- 
native approach to performing the computation would be to use an automatic differentiation 
package to evaluate numerically the required derivatives. Such an approach would mean that 
the complete algorithm could be implemented in, say, Fortran 77 or Fortran 90, thus avoiding 
both the long symbolic calculations reported on in section 6, and also the need to interface 
the Fortran code to the Mathematica package. We intend to address this issue in a future 
publication. 

( f ) Provably correct computations There is a considerable interest in performing provably 
correct computations. We note that the symbolic algorithm described in section 6 provides, 
not only an estimate of the L'^[0, oo) solution at some point A > 0, but also a provably correct 
bound on the error at A. This information could be used as input data to an interval-based 
ordinary differential equation solver, and the computation of the M(A) spectral matrix could 
be performed in a provably correct manner, giving precise information on the numerical errors 
involved. The effectiveness of this approach will need further investigation. 
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